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We study a continuum model of the lipid bilayer based on minimizing 
the free energy of a mixtu re of water and lipid molec ules. This paper 
extends previous work bv lBlom and Peletier ( 2004f ) in the following 



ways, (a) It formulates a more physical model of the hydrophobic ef- 
fect to facilitate connections with microscale simulations, (b) It clar- 
ifies the meaning of the model parameters, (c) It outlines a method 
for determining parameter values so that physically-realistic bilayer 
density profiles can be obtained, for example for use in macroscale 
simulations. Points (a)-(c) suggest that the model has potential to 
robustly connect some micro- and macroscale levels of multiscale 
blood flow simulations. The mathematical modelling in point (a) 
is based upon a consideration of the underlying physics of inter- 
molecular forces. The governing equations thus obtained are mini- 
mized by gradient flows via a novel numerical approach; this enables 
point (b). The numerical results are shown to behave physically 
in terms of the effect of background concentration, in contrast to 
the earlier model which is shown here to not display the expected 
behaviour. A "short-tail" approximation of the lipid molecules also 
gives an analytical tool which yields critical values of some parame- 
ters under certain conditions. Point (c) involves the first quantitative 
comparison of the numerical data with physical experimental results. 



1 Introduction 



We s tudy the potential utility of a continuum paradigm (jBlom and Peletier 



20041) of the lipid bilayer as a mesoscale filter in multiscale simulations of blood 



flow. The two main points are first that we introduce a more physical model of 
the hydrophobic effect into the paradigm (for consistency with microscale simu- 
lations), and second that we investigate the role and meaning of the paradigm's 
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parameters (to connect with macroscale simulations) by performing the first 
quantitative comparison of numerical solutions of the paradigm with physical 
experimental data, and in so doing provide a method for determining parameter 
values. 

As is well-known, the cell is the fundamental element of all living matter. The 
activity of the cell sustains life and the cell itself is sustained by a metabolism 
which utilizes the mass transfer through its membrane. The cell membrane is 
composed of lipid molecules with proteins and other components floating in it. 
The dynamics of this lipid bilayer membrane become especially important in 
the case of dispersed components in the blood, such as red blood cells, white 
blood cells, platelets and so on, because the deformation dynamics of these 
membranes directly affect the mass transfer in the blood. These membranes 
are often modelled as hyperelastic due to the presence of a cytoskeleton. On 
the other hand, a liposome — which is composed of lipid bilayers only — is 
modelled as a two-dimensional fluid membrane, because the membrane lipid 
molecules can easily move laterally within the bilayer. Liposomes are used 
as drug delivery agents and artificial oxygen-carriers in blood. Nevertheless, 
all these cases are defined by the lipid bilayer and the behavior of the bilayer 
is responsible for the mass transfer through the cell membranes. Hence, the 
modelling of the lipid bilayer membrane from a molecular level through to the 
contin uum level is highly anticipated to predict the mass transfer behavior in 
blood jMchcdlishvili and Maedal. l200lb iBorvczko et~ail l2003t ISugii et all 120051: 



Vaidva et aUl2007t) . 



Blood flow is a complex and major aspect of the growing field of mu l tiscale 
simulations of the hum an body ( Avton et al. . 20021 : Sugii et all . 120071 : iBoall . 



20021 iMouritseil l2005h . It is multiphase, with red blood cells (RBCs), white 



blood cells, platelets, artificial drug delivery agents (DDAs) and other compo- 
nents dispersed in a water-like plasma. It involves multiphysics, with blood 
components having viscoelastic properties and mutually interacting while also 
exchanging mass and undergoing chemical reactions. Moreover, it is multiscale; 
quantum mechanics at the smallest length and time scales upwards through 
other physics at intermediate scales, and continuum mechanics at the longest 
scales are all required to properly understand, predict, and control the behaviour 
and properties of blood. Advances in computer processing speed and storage 
capacity have made multiscale simulations possible, highlighting the need to 
understand the computational and physical interactions between the scales. 

Here we focus on one crucial intermediate scale, or mesoscale, at which 
both molecular physics and continuum mechanics are important. This is at the 
membrane which surrounds RBCs and some DDAs, being typically only two 
molecules in thickness but extending laterally for several micrometres. Under- 
standing how membrane composition affects deformability, and how deformation 
affects the mass-transfer properties of RBCs and DDAs, are key to the multi- 
scale modelling of blood flow, since the blood vessels in which occur the greatest 
mass transfer from these components are also those in which they undergo the 
greatest deformations due to large shear forces and mechanical constriction in 
narrow capillaries. One significant challenge for multiscale simulations almost 
independent of large increases in computational power and storage, is the need 
to reliably transf e r info rmation from one level of the simulation to another. 
Presently ( Yound . l200ll ). this is mostly done ad hoc and offline, and may omit 
crucial mesoscale dynamics. 
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The membrane is a bilayer of lipid molecules (lipids) (jMouritsen , 20051) . 
The dipolar charge distribution in the lipid "head" groups make them able to 
bond with water molecules, while "tails" are indifferent to the presence of water. 
This amphiphilic nature leads to t he hydrophobic effect which physically forms 
the bilayer and gives it integrity (jChandlerl . 12005). In this paper, we consider 
bila yers composed of one typ e of lipid. 

iBlom and Peletierl (|2004!) base their continuum paradi gm, here i n call ed the 
"BP paradigm", on the mesoscopic dynamics framework of [Fraaiiej (|l993[) , min- 
imizing a free energy for a system of lipid and water molecules. Formally, the 
intrinsic free energy of the system is minimized with respect to a constraint 
that the (unobservable) distribution of the molecules generates the (observable) 
continuous volume fractions, thus assuming that the microstate has relaxed to 
equilibrium over the relatively long time scale of the continuous description. 

We introduce a new model of the inter-molecular interactions into the BP 
paradigm which in contrast to the original model closely represents the un- 
derlying physics of the hydrophobic effect. A further improvement is a term 
(3 to control the decay rate of the interaction strengths. Our new model and 
other improvements anticipate being able to determine the values of more of 
the model parameters from microscale simulations (or first principles). Since 
tests not reported here show that both models have very similar computational 
costs, we concentrate largely on our new model. Our point is not that our model 
has computational advantages or performs better under certain conditions, but 
rather that it has a more physical basis and so enables a better connection to 
microscale simulations or first principles. It is therefore encouraging that our 
more physical model give s results qualitatively similar to the original model of 



Blom and Peletierl (|200j). Moreover, we show here that a key characteristic of 



the bilayers produced by the new model depends on a model parameter in a 
physically realistic way, in contrast to the original model ( H3.2.ip . 

Understanding the paradigm's parameters and the numerical roles they play 
is important for connecting to macroscale simulations and physical experiments. 
By considering the analytical origins of the paradigm, we show that some of 
those parameters which are at first sight physical (as opposed to purely numer- 
ical) are seen to be largely numerical, in that they cannot be directly connected 
a priori with physical measurements. Properties of numerical bilayers must be 
compared a posteriori with physical properties in order to set some parameter 
values. The numerical bilayers are obtained by a method of solution new to 
the paradigm, and are for the first time quantitatively compared with physical 
bilayers. 

In more detail, the system of water molecules and heads and tails of lipids has 
a free energy split into an ideal part roughly corresponding to the Helmholtz free 
energy, involving only connectivity interactions, and a non-ideal part represent- 
ing inter-molecular int eractions. Lipid structure and configuration are therefore 
explicitly represented. IBlom and Peletierl ( 20041 ) formulate a non-ideal part of 
the free energy with a term reflecting the (global) compressibility of the sys- 
tem and another modelling the (local) hydrophobic interactions; the modelling 
of the hydrophobic interact i on ter m is an open question and is not inherent 
to the formalism of iFraaiiel (|l993f ). In this paper we introduce a new model 
of the hydrophobic interactions which captures the physics underlying the hy- 
drophobic effect responsible for the formation and integrity of ordered states of 
amphiphiles, based on the following discussion. 
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Liquid water is a dynamic hydrogen bond network in which each water 
molecule forms up to four hydrogen bonds with its neighbours. The non-zero 
dipole moment of lipid head groups makes them able to accept hydrogen bonds 
from water molecules (but unable to donate a bond to each other): they are hy- 
drophilic. By contrast, the hydrophobic lipid tail groups are unable to form hy- 
drogen bonds, although thermodynamic and electrostatic interactions between 
water molecules and tail groups occur, but in liquid water at room temperature 
the hydrogen -bond energy is typically an order of magnitude stronger than such 
interactions ( Immergut . 1991 : MouritsenL 20051 ). 

A cavity with a structured "surface" in the hydrogen bond network forms 
around a hydrophobic moeity, causing a decrease in the entropy of the system 
( Kronberg et al. . 19951 ). An entropic force acts to gather together hydrophobic 
moieties so as to minimize the disruption to the hydrogen bond network. The 
physical origin of this hydrophobic effect is that water molecules close to a 
sufficiently large hydrophobic moiety no longer participate in four hydrogen 
bonds; with no attractive force towards the hydrophobic moeity, these molecules' 
remaining bonds now draw them away from the moiety. It is thus because lipid 
head groups can be nodes in the hydrogen bond network while tail groups cannot 
that bilayers and other structured lipid assemblies form, and this is the basis of 
our model. 

Of the physical parameters, we take the system to be incompressible, leaving 
the effects of the compressibility parameter p to future work. The relative 
tendency of heads and waters to form hydrogen bonds is modelled here by 
the new parameter 7 which is set as unity in this paper: the effects of 7 are 
also left to future work, while we simply note here that 7 moves the paradigm 
beyond modelling heads as attached water molecules. Herein we investigate 
the effects of the key physical parameters a, e, /3, Co, and m. Respectively these 
represent temperature effects, lipid head-tail group separation distance, decay of 
the interaction strength, and a "background co ncentration" and "excess m ass" 
of lipids (these last two terms, inherited from Blom and Peletier (|2004 ). are 
clarified in this paper). 

The paper is structured as follows. We introduce into the BP paradigm a 
new model of the inter-molecular interactions in ij2.ll which forms the basis of 
the numerical solutions. The parameters and lipid model are discussed here. 
Euler-Lagrange equations, whose solution minimizes the free energy functional, 
are derived in §2.2) and a novel numerical approach to solving them is given 
in £ 12. 3[ along with a sample numerical result. All solutions and discussions are 
based on a one-dimensional model. The "smoothed" nature of the paradigm and 
the choice of lipid model are discussed in ^2.41 

§T|] connects the BP paradigm to physical in vitro measurements of lipid bi- 
layers. The summary of bilayer properties in . II is used in §3.2l to describe how 
numerical solutions can be calibrated to physical data. Numerical solutions are 
calibrated in this way for the parameters of interest in § §3.2.11 2.3. In particu- 
lar, a short-range interaction (or "short-tail") approximation based on the new 
parameter (3 is introduced in ^3.2. 21 and its effects studied analytically and nu- 
merically. This work enables guidelines on the choice of parameter values to be 
given in in §3.31 The conlusions are in §3 
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2 The new model: derivation and numerics 
2.1 The new model 

The original model of the hydrophobic interaction acted to move tails away from 
heads and waters by penalising proximity between them, mimicking the effect 
of the hydrophobic force but not the underlying cause, namely the uneven force 
distribution on water molecules in close proximity to lipid tails described above. 
Our approach, in direct contrast to the original model, is to promote water- water 
and water-head (but not head-head) proximity, modelling the hydrogen bond 
network, and effectively ignore the hydrophobic tails in direct imitation of the 
underlying physics. The relative strength of the water-water bonding preference 
to the water-head bonding preference is controlled by a par ameter 7, meaning 



that h eads are no longer the unphysical "attached waters" of iBlom and Peletier 



(120041) 



Formally, our system comprises "waters", each represented by a single "bead", 
and "lipids", each represented by a "head" bead and a "tail" beacfl connected by 
a rigid massless rod of length e. The one-dimensional model has two lipid groups 
aligned in the ^-direction, normal to the bilayer plane; one group whose tails, 
having normalized density u(x), point in the positive x-direction, and the other 
with tails of normalized density v{x) pointing in the negative cc-direction. The 
head beads of the first group have normalized density r_ e w(x) = u(x + e), and 
similarly for the second. Water beads have normalized density w(x). The lipids 
have here been chosen as the simplest allowed in the paradigm; we later ( §2.4j) 
argue that more complex models mostly would not improve the utility of the 
model as a mesoscale "filter" in a multiscale simulation. 

The total free energy of the system consists of three parts in the form 

E = T [t](u) + r](v) + r)(w)] + — (1 — u — v — T~ e u — T e v — w) 2 

f {2A) 
+a / wk * [w + 7(r_ e w + r e v)] . 



The new model differs from the original in the third part; the meanings of all 
three parts, the parameters and variables follow. 

The first term, favouring spreading, represents the entropy of the system, in 
which rj(s) = slogs for non-negative s and r](s) = 00 otherwise, and where T 
is the temperature of the system. The second term is a potential energy due to 
compressibility, where p is the system pressure. 

The new form of the third term, modelling directly the underlying physics 
of the hydrophobic effect, involves a water- water term wk * w and a water- head 
term wk * 7(r_ e w + t € v), where * indicates convolution in the form 

(fk * g){x) = { f(x)k(x - y)g(y)dy . (2.2) 



The overall strength of these interactions is controlled by a with their relative 
strength controlled by 7. The interaction kernel takes the form 

k(s) = kq — k(s) for k(s) = Sp(s) (2-3) 



1 A11 beads are of zero dimension. 
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where <5^(s) is a general smooth function with the properties 

^(±00) = 0, fs p = l, (2.4) 



and the constant kq is chosen so that J k = 1. In this paper, we define k as 

<») = W = ^ 1 (2-5) 

although other choices could be considered. The new kernel k(s) rewards prox- 
imity and thus represents an a ttractive water- wat e r and water- head force, in con- 
trast to the original model of Blom and Peletier (2004) which penalised water- 



tail and head-tail proximity. The new parameter (3 controls the decay of the hy- 
drophobic interaction, and will be shown in §3.2.2l to introduce a straightforward 
a nalytical tool. (The notati on was chosen so that the original interaction kernel 



of iBlom and Peletierl (|2004h is a special case of ours, namely k(s) — kq — k(s) 
when (9 = 1, but the hydrophobic integrand is quite different.) 

The system is taken to be infinite with an averaged de nsity cn for both u 



and v . The original question addressed to the paradigm in blom and Peletier 



(|2004l ) was, in effect, whether the minimizer of the energy functional E favors 
aggregation or spreading when an excess mass m of lipids is added to the system, 
with all other parameters fixed. Here, we concentrate on the utility of the BP 
paradigm, showing that physically realistic bilayer density profiles can be ob- 
tained from it numerically. We mainly focus on the results from our new model, 
for the reasons given in the introduction. To simplify the analysis and compu- 
tation we consider the case of periodic cells of length 2L. All the discussions 
are valid for the infinite system. 

2.2 Derivation of the Euler-Lagrange equations 

Here we take the system to be incompressible, p — 00, so that 

1 — u — v — T- t u — t c v — w = . (2-6) 
The energy functional can be simplified as 

[t)(u) + T)(v)] + ct f (1 — u — v — T- e u — T e v)k * (1 — u — v) (2.7) 

subject to constraints 

1 — u — v — T- t u — r e v — w ^ , (2.8a) 

(u + v-2c )=m. (2.8b) 



The result of scaling T into a is that we can consider the temperature effects by 
varying a (i j3.2.3[) . We have also dropped the entropy of the water molecules, 
which is justified since although the entropy changes of the water associated with 
reduced configurational arrangements ar ound hydrophob i c moie ties actually as- 
sists solvation, the effect is very small (jKronberg et al. . 1995 ). Furthermore, 



we have taken 7=1, effectively indicating that the electronegativities of wa- 
ters and heads are equal, and leaving the effects of the relative strength of the 
water-water to water-head bonding preference to future work. 
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Using the definition (|2.5[) of the interaction kernel and the mass conservation 
constraint (|2.8b[) the energy functional becomes 



Er 



m \ 



(1 — U — V — T^ e U — T e v) 



[r)(u)+r)(v)]+a(l-2c - . }/ , 

—a / (1 — u — v — T- e u — t € v)k * (1 — u — v) 



(2.9) 



where we have chosen kq — (2 — e L ^)/2L. Using the method of Lagrange 
multipliers we rewrite the energy functional as 



Et — Ej 



K 



M 2 + h 



m — / (u + v — 2cq) 



(u + v — 2c ) 



(2.10) 

where K and A± are Lagrange multipliers and /i = (u + v + T- e u + r e v — 1)+, 
with (-) + = max{-,0}. 

Carrying out calculus of variations in a formal way, assuming that the order 
of integrations and translations can be changed wherever necessary, we derive 
the Euler-Lagrange equations 



= logu — an * (2u + 2v + 2r_ e w + T_ e w + r t v) + Kji + Kji(x + e) + A , 



(2.11a) 

= logw — an * (2u + 2v + T- C u + t € u + 2r e v) + K[i + K[i{x — e) + A , 

(2.11b) 

where 

A = A_ - A+ + 1 + 3a - 2a(l - 2c - m/2L) . (2.12) 

We solve the Euler-Lagrange equations by first replacing them with evolution 
equations based on gradient flows: 

u t = — log u + an * (2u + 2v + 2t_ £ u + r_ e v + T e v) — K[i — K[i(x + e) — A , 

(2.13a) 

v t = — log v + oik * (2u + 2v + T_ e ii + T f u + 2t c v) — K[i — Kfjb[x — e) — A , 

(2.13b) 

for the gradients Ut = —(6E/Su), v t = —(8E/5v), with given values of K, A. 
These equations are solved numerically in the next section. 

There are thus nine model parameters in total: seven apparently physical 
and two (K, A) strictly numerical. K and A will be chosen subject to stabil- 
ity, symmetry, and water conservation considerations described shortly. Of the 
physical parameters we consider herein e, (3, Co, m, a, as described previously. 

2.3 Numerical scheme and solutions 

Our numerical scheme solves the gradient flow equations (|2.13al b) to find the 
densities u, v which minimize the Euler-Lagrange equations (|2.11al b). The ap- 
proach uses a discrctizcd grid with finite difference formulae for Ut,Vt- Both 
first- and second-order backward differencing was used, with the results in 
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Figure 1: A sample "bilayer" for the parameter set a = 3,e = 2,(3 = l,co = 
0.024, to = 0.05 * 2L. 



close agreement. Generally, very small time steps were required. The nu- 
merical domain w a s tak en to be periodic and of length 2L, in contrast to 
Blom and Peletierl ( 20041 ) who used a finite domain with small decay at the 



edges. In both approaches, u, v can and do deviate from Co at the boundaries. 
In tests using the orig i nal m odel, our numerical scheme reproduced the profiles 
of iBlom and Peletierl (j2004T ). as near as we can tell, given that their simula- 
tion conditions were not fully specified. Typical runtime on a machine with 
two 2GHz AMD Opteron 270 Dual Core procesors with 16GB of DDRS-667 
SDRAM is four minutes. 

The values of the Lagrange multipliers K, A are not specified by the model. 
Indeed, we require an explicit penalty term in the algorithm, effectively replacing 
A_ — A + in A with X*(J(u + v — 2c ) — m). Within the range of values of K, A* 
for which the numerics are stable, K is chosen large enough to ensure w 
but no larger, and A* is chosen to ensure that the solutions are symmetrical, as 
expected. Typically, K is of the order of 10 3 (but can be as large as O(10 4 )), 
whereas A* is around 10 2 . 

A sample result is shown in figure [1] The system has separated into a well- 
defined bilayer-like profile, in which the hydrophobic tail region is separated 
from the water region by two peaks in the hydrophilic head group density. 
Although it is not our purpose to compare solutions of our model with those 
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of the original model of iBlom and Peletierl (|2004l ). we here note briefly that the 



forms are similar, but with lipids being slightly more strongly drawn into the 
bilayer of the new model. All the results presented in this paper are independent 
of domain length 2L (above a certain value) and grid resolution (beyond some 
level of coarseness). 

2.4 Smoothing and the choice of lipid model 

Referring to our discussion in §T]of the formulation of the BP paradigm, the den- 
sities of physical lipid moieties are smoothed spatially and temporally into the 
densities of the model "lipids". In general, there will not be a one-to-one corre- 
spondence between each such "lipid" and the physical lipid molecules, but there 
is a rigorous correspondence between the densities of the one and those of the 
other. We prefer to call the paradigm's "lipid" a model lipid component. These 
components represent an as-yet unclear spatial smoothing (coarse-graining at 
the same length scale) and temporal smoothing of the molecular-level informa- 
tion. Regardless, we are free to choose any lipid component model; we have here 
picked the simplest. This does not mean that any information from the physical 
lipid is lost, but that it is smoothed into (some sum of the densities of) the two 
beads per component which we have. This is important from the perspective of 
integrating this paradigm into a multiscale simulation, for which we are inter- 
ested in general bilayer characteristics such as bending rigidity which rely more 
on physical thickness rather than details of the bilayer profile. A simple rod- 
and-two-beads model does the job of capturing all of the molecular data into a 
much smaller number of variables, and in general no more complex model lipid 
components need be considered. One exception might be to replace the rigid 
rod by a spring, since especially in the pr esen c e of e mbedded proteins physical 
lipid molecules can stretch considerably ( Led . 2003 ). However, even this case 



may be addressed by the current model lipid component. 

The temporal smoothing occurs in moving from the discrete to the contin- 
uous description, working at the lipid length scale. The continuum densities 
are identified formally with the discrete ones over a sufficiently long time scale, 
but it is not clear what this time scale is. This smoothing, although unclear in 
detail, is central to the smoothed paradigm and is important not just because it 
enables access to longer time scales and hence can act as a bridge from a lower 
level simulation to a higher level one, but also because it captures some detail 
of the real thermal motion of the lipid molecules within the bilayer, which is a 
defining characteristic of bilayers and plays a key role in their function. 

Together, the spatial and temporal smoothings are evident in the numerical 
results in that the model lipid components have combined to create a bilayer of 
total width greater than the naively-expected 2e, as can be seen in figure Q] (see 
a more formal discussion in 



3 Connecting the paradigm with the physics 
3.1 Physical bilayer properties 

Three biologically significant membrane characteristics are (1) the el astic mod- 
uli, (2) the intrinsic monolayer curvature, and (3) the bilayer thickness ( MouritsenL 
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Figure 2: The e xperimentally-determined structure of a DOPC bilayer (figure 
reproduced from I Wiener and White! (Il992l ) with kind permission of the authors 
and the Biophysical Society.) 



20051 ). The latter is the focus of the present section ((1) and (2) require working 



in higher dimensions). 

The averaged thickness g?sz of the hydrophobic core, or saturation zone, is a 
common physical measure of bilay er thickness. In a physical system this can be 
increased by the following means (lMouritsenl . [2005l . §8.3): increasing the length 
Irof the tails; replacing the double carbon bonds by single bonds in the tails; 
decreasing the degree of hydration; increasing the cholesterol concentration; 
decreasing the temperature. (The latter means each lipid has less kinetic energy, 
making all lipids on average less convoluted, and hence longer on average, while 
in a bulk effect lowering the temperature makes the arrangement of lipids more 
crystalli ne. Both effects in crease dsz-) These effects can increase dsz by several 
percent (|Mouritsenl l2005l §9.2). 

The main approaches used in physical experiments to determine the time- 
averaged structure of lipid bilayers, and hence their thickness, exploit the high 
structural periodicity in the direction normal to the bilayer, for example in com- 
binin g diffraction data from x-ray and neutron scattering ([Wiener and White! . 



19921 ). One such data set is represented in figure[5]for the DOPC lipid molecule. 



In this figure, the water density measures only the waters of hydration (those 
bonded to head groups). 
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3.2 Making comparisons with physical meaurements 

For our analysis the key features of figure [5] are the saturation zone, defined as 
that region where w « 0, the density curves of the end of the tail (here, the CH3 
moiety of DOPC), and those of the head group (here, choline and phosphate). 
The DOPC bilayer will form the basis of the remaining work in this paper, 
but the key point is that any single-species bilayer can be simulated by the BP 
paradigm, when basic structural details are known from experiments. 

We consider only those numerical results with a clear bilayer structure like 
that shown in figure [TJ namely with a single saturation zone of width dsz, a 
single tail peak with exactly two transitions from concavity to convexity, and two 
head peaks of equal width c?hz , each likewise with exactly two transitions from 
concavity to convexity. Macro-separation of the two species u and v can and 
does occur without these criteria being met, but our analysis is here concerned 
only with those profiles rigourously of this form. The head zone width dwz is 
defined in the following way. Using the left-hand head peak, let xl, xr be points 
immediately to the left and right of the peak satisfying h x — 0,h xx > 0, i.e. 
local minima, with the local maxima of the head peak located at x* . We then 
findxi,a;2 from h{xi^) = \{h{xL,R) + h{x*)) and define g?hz = x^—xy. Defining 
the head zone width in this way as running from th e left-hand midpo int to the 
right-hand midpoint of the peak is not unusual (e.g. Mouritsenl ( 2005 . fig 8.1)), 



but none of our conclusions is changed significantly by defining it as xr — xl- 

Lipid bilayers have different cZhz and dsz depending on the lipids of which 
they are composed. Specifically, the ratio (fez /dsz depends on the choice of lipid 
molecule (other factors such as temperature being equal) and so characterizes 
the bilayer properties for our purposes. For the DOPC bilayer of figure [2j 
daz/dsz ~ 0.4. 

We now show that varying the key paradigm parameters of i j2.ll enables a 
solution to be found corresponding closely to any desired physical bilayer, with 
the DOPC bilayer as our example. The method of selecting the values of the 
parameters is summarized in 



3.2.1 Co, m and the critical micelle concentration 



Blom and Peletierl (|2004r ) called Co the "background concentration" and m the 
"excess lipids" above the background concentration. This terminology stems 
from the physical phenomenon known as the critical micelle concentration (CMC), 
in which the monomer density of lipids in solution only increases up to the CMC, 
beyond which the excess lipids aggregate into o rdered stru ctures. The type of 



ordered structure depends on the lipid geometry (|Boall . l2002l ). The BP paradigm 



takes Co to be the background concentration to which an excess quantity of lipids 
m is added. 

When implementing the paradigm numerically we face the problem that the 
smoothing of §2.41 only reveals the lipid species a posteriori once the numerical 
profiles are compared with experimental ones. Without knowing the species, 
and without also knowing a priori the temperature of the system, the CMC 
cannot be known at the start of the simulation, and hence neither can Co and 
m. Moreover, the CMC refers to an average lipid density, which here would be 
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Figure 3: duz and dsz against Co for three different pairs (ce,e), with all other 
parameters equal. Straight lines of best fit have been drawn. 



Consequently, the true background density is m/L + 4cq, where Cq is the value 
of Co for which aggregates first form (with all other parameters fixed), and the 
total excess of lipids is 8L(co — Cq). Since m and Co combine to form the CMC, 
we can numerically fix m at a working value and vary Co only. 

We note here that the above argument suggests that u and v should decay 
to (m/L + 4cq)/4 far from the bilayer, but both the periodic domain used in 
the present numerical si mulations, and the finite arbitrary value of L used in 
the numerical method of iBlom and Peletierl (|2004l ). allow u and v to vary away 
from the background concentration at the "edges". 

In figures [3] and 0] we consider three data sets in which the parameters a,e, 
and Co vary, the values of a and e being given in parentheses on the graphs. 
The widths g?hz and dgz are normalized on e in these figures. Numerical data 
is represented by points and a straight line of best fit is drawn in each case. 

In the numerical system, as Co increases beyond Cq, the excess lipids thus 
supplied should be drawn into the bilayer with the numerical background con- 
centration remaining more or less the same. This required increase of c£hz and 
dgz with increasing Co can be seen in figure [3] Although this is encouraging, 
this one-dimensional study cannot yet reveal the true physicality of the mod- 
els in terms of varying Co. In a physical system, as more lipids are added to 
the solution beyond the CMC, more micelles (or othe r ordered st ructures) are 
formed, with the average size of each barely affected (jBoal 2002). Such a test 
must wait for two-dimensional results. 

More importantly, the data can be combined into the ratio dnz/dsz as 
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Figure 4: The ratio c^hz/^sz against cq for the three parameter sets of figure [3] 
Straight lines of best fit have been drawn. 



shown in figure HJ Taking the example of a DOPC bilayer, the data shows 
that we can choose suitable parameter sets such that c?hz/^sz ~ 0.4, namely 
(a,e,co) = (3,2,0.028) for which d H z/d S z = 0.41, (a,e,c ) = (4,2,0.022) for 
which dnz/dsz = 0.37, and (a, e, c ) = (3,3,0.052) for which dnz/dsz — 0.42. 
We expect that ultimately a can be set physically (see also ^3.2.30 and that the 
role of e is more one of clarifying the structure 0j3.2.2j) . so that effectively here 
our only choice would be Co, and the data shows that a value can be chosen 
which yields a bilayer characteristic close to the desired one. This is expected 
to hold for lipids other than DOPC since there was nothing special about our 
choice of this molecule. 

Although it is not our intention to directly compare the new model with 
the original, it is worth noting here that the original model does not show the 
expected increase of g?hz with cq: indeed, an inverse relationship between c?hz 
and Co holds for numerical data generated by the original model, as shown in 
figure 



3.2.2 e and the interaction decay length 

That e is not an actual lipid length has already been discussed. What role, then, 
does it play? 

Varying e changes the degree of separation of the head and tail regions, larger 
e giving clearer bilayer structure. With our definition fp. llll) of what constitutes 
a realistic bilayer profile, we need only consider a few order unity values of e; 
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04 0.045 0.05 0.055 



Figure 5: Example data for the original model: g?hz against cq for (a, e) = (4, 2), 
plotted on the same scale as figure [3] A straight line of best fit has been drawn, 
indicating a non-physical inverse relationship between cZhz and cq. 



see JO] 

Turning to the parameter /3, which controls the interaction decay length, 
we consider another special case of the kernel function of equation (|2.5|) . Our 
"short-range interaction" approximation takes (3 — ► 0, so that the kernel function 
approaches the delta function 



k(s) 



s = ; 



(3.2) 



The Euler-Lagrange equations (|2.11|) then become 



= log u - a(2u + 2v + 2r_ £ u + r^ e v + r e v) + K[i + K[i{x + e) + A , (3.3a) 
= log v - a(2u + 2v + T- £ u + r e u + 2r e u) + Kfj, + Kfi{x - e) + A , (3.3b) 



which in the case of K — > oo reduce to 



logw — a(2u + 2v + 2r_ e u + r e v + r e v) = —A , 
logv — a(2u + 2v + r_ e u + r e u + 2r e v) = —A . 



(3.4a) 
(3.4b) 



Separating A into a constant term plus a term dependent on c by rewriting 
(|2.12p as A = A + 2a(2c + m/2L), we look for constant solutions u = v = c, 
obtaining 

logc — 4ac = — A. (3.5) 
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Differentiating (|3. 5|) with respect to c yields the critical concentration c c as the 
first necessary condition for the existence of a solution: 

c = c c = ^- . (3.6) 

The sho rt-range interaction tool t hus gives a simpler way to find the same result 
(|3.6p as lBlom and Peletier ([2004) . To simplify our discussion, we assume that 
Co = c c — 1 /4a and take a ^ 1. 

The other necessary condition for the existence of a solution is estimated 
as follows. From the incompressibility condition (|2.8a|) we see that c ^ 0.5. 
Substituting c = 0.5 into (|3.5[) gives 

A min = -log 2 -2a . (3.7) 



When there is phase separation, (13. 5|) must have more than one solution: indeed, 
A must also be less than A max where 

Amax = logc c - 4ac c = - log(4a) - 1 , (3.8) 
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in which case there is only one solution. When (3 — the incompressibility 
condition implies c ^ 0.25, and therefore 



A min = - log 4 



(3.9) 



With A m i n < A < A max there are many possible states. The relevant question is 
whether there exists a global minimum of the <5-function free energy 



u + v — 2cq — —J (1 



- U — V — T_ e U - 



r e v) 



E D = lim Ei -- 

(3-10) 

With A chosen between its minimum and maximum values, we assume that 
u = v = c\ for —l<x<l and u = v = c 2 elsewhere, where c\ and c 2 are two 
distinct solutions of (|3.5p . We also assume that the solution is periodic with 
period 2L. The mass conservation (|2.8bj) gives 



+ 4L(c - c 2 ) 
4(ci - c 2 ) 



(3.11) 



which in p.lOp yields 



Ed 
4 



(7YI \ 
2(ci + c 2 ) + 2c + —J (ci - c 2 ) 



+L 



c 2 logc 2 + a (c + ^\ (4c 2 - 1) + 4ac2 - ea(ci - c 2 ) 2 



(3.12) 



In figure O we have plotted Ejy as a function of A, for (a,e) — (4,2). We can 
see that the global minimum occurs when A = A max . 

Returning to the question of finding a desired profile, figure [7| shows how 
^HZj^sz and their ratio vary with (3 for all other parameters fixed. A multi- 
lammcllar profile began to appear for [3 < 0.75, while separation without a 
saturation zone appeared for (3 > 2.25. Crucially, a bilayer-like profile could be 
found for all values of (3 in the given range, meaning that the important ratio 
d-Hz/dsz can be fine-tuned to the desired accuracy. See 



3.2.3 a and the effects of temperature 

Recalling the physical factors influencing dsz listed in §3.1[ since a is inversely 
proportional to T, increasing a should increase the measurable dsz- Because 
this can be seen in figure^ it appears that temperature-related mechanisms can 
be captured in ID. Indeed, turning to the incompressible free energy functional 
(|2.7p in which a is scaled on T, increasing a (decreasing T) reduces the effect of 
the entropy relative to the interaction terms. This makes physical sense in that 
with less kinetic energy the hydrogen-bond network is more strongly preserved. 

Further, the other physical temperature-related mechanism is that by which 
lower temperatures straighten the lipid tails on average. (Although the head 
groups also become "straighter" in the sense that they change position less and 
therefore have a longer averaged shape, the main effect is in the longer and more 
deformable tails.) As a result, we would expect the ratio d^z/dsz of our model 
to decrease as a increases (T decreases), and this is seen in figure 2J 

Finally, we caution against one potential method of setting a. That increas- 
ing a (decreasing T) thickens the physical bilayer (in terms of c?sz) is from the 
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Figure 7: c?hz, dsz, cfez/dsz for different values of j3. The other parameters 
were fixed: in particular, (a,e, Co) = (3,2,0.024). For c?hz the drawn line of 
best fit is straight ax + c, for dsz it is of the form ax + b/x + c, and for the ratio 
it is exponential and of the form aexp(6x) + c. 



view point of moving with the bilayer as it undergoes thermal fluctuations. But 
increasing a (decreasing T) should suppress these thermal fluctuations. If the 
model were capturing these undulations as well then we might expect dsz to 
decrease as a increased (T decreased) because the order of magnitude of the 
thermal fluctuations is larger than that of the changes in dsz (a few percent) . 
Since we see the opposite, this precludes using the amplitude of an averaged 
thermal undulation to fix a. 

3.3 Choosing a parameter set 

Based on §£ !3.2.1|3.2.2|3.2.3t we can set the parameters to reproduce the pro- 
file of a desired lipid species bilayer whose properties are known a priori from 
experimental data as follows. 

Choose a = 3 or 4, noting from figure |J] that from the numerical viewpoint 
a controls the sensitivity of the ratio dnz/^sz to variation in cq. Then, recalling 
that the effects of e and (3 combine, we pick e = 2 or 3, essentially only requiring 
clear separation of the head and tail regions. Then we set /3 = 1 to start. Next, 
since Co , m combine with m fixed and cq to vary, we run the numerics with a few 
Co until G?Hz/cfez is close to the desired value. These results are then fine-tuned 
by varying /3, which smoothly varies the values g?hz and dsz- 
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4 Conclusions 

We introduced a more physical mode l of the hydrophobic effect into the contin- 
uum paradigm of iBlom and Peletier and showed that one-dimensional 



numerical solutions can reproduce key characteristics of physical lipid bilay- 
ers. The realised goal was to show that a more physical model enabling better 
connections with microscale numerics and first principles is capable within the 
paradigm of generating physically realistic bilayer density profiles. In particu- 
lar, the mechanically-important bilayer thickness can be reproduced. Indeed, 
various key characteristics were shown to behave physically, at least one in con- 
trast to the original model. The coarse-graining or smoothing of the molecular 
information inherent in the paradigm was clarified, and the smoothed numer- 
ical data calibrated to measurements of physical bilayers. Examining some of 
the key parameters in turn, we gave a strategy for setting their values, noting 
that future work, especially in higher dimensions, could make this process even 
more robust by further appeal to physical arguments. In particular, a is for- 
mally linked to (the inverse of) the temperature, and has been shown here to 
have that effect on the numerical results, while f3 should be related to the range 
of the hydrogen bonding forces. Already, however, the new parameter j3 has 
been seen here to introduce the short-range interaction tool, allowing analytical 
results to be obtained with greater ease. 

Much work remains to show that the BP paradigm is a viable mesoscale 
tool for multiscale simulations. The main aim is to work in higher dimensions, 
and include compressibility effects by allowing p to vary. The former represents 
a significant computational challenge. If the conclusions of the current paper 
are supported by higher-dimensional work, then the paradigm can be made 
as physical as desired by varying 7, considering other k (and indeed whether 
"promoting proximity" is the best model of the underlying physics), including 
electrostatic effects, weak head-head van der Waals repulsion terms, and so on. 
The level of detail will rest on computational issues, in particular cost-benefit 
considerations in the light of the paradigm's use as a mesoscale numerical filter. 



Acknowledgements 

P.L.W. gratefully acknowledges the support of the Japan Society for the Promo- 
tion of Science and the University of Tokyo's Intelligent Modeling Laboratory 
during part of this work. His thanks for many helpful discussions go to mem- 
bers of the IKEMEN group of the Fluid Engineering Laboratory of the Univer- 
sity of Tokyo. We also thank the RIKEN institute for their support through 
the project Research Program for Computational Science, R&D Group for the 
Next-generation Integrated Simulation of Living Matters, which is supported by 
the Japanese Ministry of Education, Culture, Sports, Science and Technology. 



References 

Ayton, C, Smondryev, A. M., Bardenhagen, S. C, McMurty, P., and Voth, 
G. A. (2002). Calculating the bulk modulus for a lipid bilayer with nonequi- 
librium molecular dynamics simulation. Biophys. J., 82:1226-1238. 



REFERENCES 



19 



Blom, J. G. and Peletier, M. A. (2004). A continuum model of lipid bilayers. 
European J. Appl. Math., 15:487-508. 

Boal, D. (2002). Mechanics of the Cell. CUP, Cambridge. 

Boryczko, K., Dzwinel, W., and Yuen, D. A. (2003). Dynamical clustering of 
red blood cells in capillary vessels. J. Mol. Model., 9:16-33. 

Chandler, D. (2005). Interfaces and the driving force of hydrophobic assembly. 
Nature, 437:640-647. 

Fraaije, J. G. E. M. (1993). Dynamic density functional theory for microphase 
separation kinetics of block coploymer melts. J. Chem. Phys., 99:9202-9212. 

Immergut, E. H., editor (1991). Encyclopedia of Applied Physics, volume 18. 
VCH, Berlin. 

Kronberg, B., Costas, M., and Silveston, R. (1995). Thermodynamics of the 
hydrophobic effect in surfactant solutions — micellization and adsorption. 
Pure & Appl. Chem., 67:897-902. 

Lee, A. G. (2003). Lipid-protein interactions in biological membranes: a struc- 
tural perspective. Biochim. Biophys. Acta, 1612:1-40. 

Mchedlishvili, G. and Maeda, N. (2001). Blood flow structure related to red cell 
flow: determination of blood fluidity in narrow microvessels. Jpn J. Physiol., 
51:19-30. 

Mouritsen, O. (2005). Life - as a Matter of Fat. Springer, Berlin. 

Sugii, T., Takagi, S., and Matsumoto, Y. (2005). A molecular-dynamics study 
of lipid bilayers: effects of the hydrocarbon chain length on permeability. J. 
Chem. Phys., 123:184714. 

Sugii, T., Takagi, S., and Matsumoto, Y. (2007). A meso-scale analysis of lipid 
bilayers with the dissipative particle dynamics method: thermally fluctuating 
interfaces. Internat. J. Numer. Methods Fluids, 54:831-840. 

Vaidya, N. K., Huang, H., and Takagi, S. (2007). Modelling HA protein- 
mediated interaction between an influenza virus and a healthy cell: pre-fusion 
membrane deformation. Math. Med. Biol, 24:251-270. 

Wiener, M. C. and White, S. (1992). Structure of a fluid dioleoylphophatidyl- 
choline bilayer determined by joint refinement of x-ray and neutron diffraction 
data. Biophys. J, 61:434-447. 

Young, D. C. (2001). Computational Chemistry: A Practical Guide for Applying 
Techniques to Real-World Problems. John Wiley & Sons, Inc., New York. 



